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A numerical scheme for ordinary differential equations having time varying 
and nonlinear coefficients based on the state transition matrix 

Robert E. Bartels 
Abstract 

A variable order method of integrating initial value ordinary differential equations that is based 
on the state transition matrix has been developed. The method has been evaluated for linear time 
variant and nonlinear systems of equations. While it is more complex than most other methods, it 
produces exact solutions at arbitrary time step size when the time variation of the system can be 
modeled exactly by a polynomial. Solutions to several nonlinear problems exhibiting chaotic 
behavior have been computed. Accuracy of the method has been demonstrated by comparison 
with an exact solution and with solutions obtained by established methods. 


Section 1.0. Introduction. Numerical methods to solve ordinary differential equations 
are at a state where efficient and accurate solutions are routine. Accuracy is often 
achieved by adapting time step size, which in some applications is inconvenient or 
impossible. It may not be possible to adapt time step size when the equations of 
dynamics are integrated while coupled to another solver (e.g. computational fluid 
dynamics (CFD ) code). In other cases an arbitrarily large constant time step size may be 
desirable such as when externally supplied data is applied to a system at constant time 
intervals. 

One approach to address this dilemma has been the construction of numerical schemes 
based on classical solutions. An exact numerical scheme based on the state transition 
matrix (STM ) is an example. One advantage of a scheme based on the STM is that it 
offers an exact solution of the system for arbitrary time step size. The ability to compute 
at arbitrary time step size is useful, for instance, when coupling a continuous system to a 
discretely sampled control system with discrete feed back. The system can be integrated 
with absolute accuracy between samples and discrete feed back with one or several steps, 
as required to capture the discrete nature of the problem. While numerical methods for 
time invariant systems based on the STM have been in use for many years 1 , only recently 
has it been applied to time variant systems. 2 The present paper proposes a discrete form 
of the STM that has application to any linear time varying (LTV) system, and some 
nonlinear time varying (NLTV) systems. 

Consider the simple first order differential equation 

x = A(t)x + B(t)u(t) (1 ) 

where x = (x l ,...,x N ) T , x n e R, Ae R NxN ,Be R ,VxA t>t 0 . Setting the dimension of 

the matrix B identical to matrix A is purely for convenience; the results are applicable to a 
matrix B of general form. In this paper the classical method of solving differential 
equations will be used to construct the numerical solution of equation ( 1 ). Between times 
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t m and t m + h, where h is time step size, a discrete approximation will be performed using 
the solution 

t m +h 

x(t m +h) = 0(t m +h,tjx(tj+ + h,T )B(T )u(T )dT (2) 

t m 

For any LTV system, the STM between times t m and t m + h is given by the Neumann 
series 


r t m +h c t m +h <• t x 

<S?(t m +h,t m ) = I + J A(t)dt + J A(t x )J A{t 2 )dt 2 dt x 
+ \‘ t m+h A(t, )j A(t 2 )j ' A(t 3 )dt 3 dt 2 dt l 4 — 


(3) 


This result is called the Peano-Baker formula. Only in special cases, such as when the 
matrix A is constant, is this solution equal to the more familiar expression 


®(t m + h, t m ) = exp 


■ t+h 


A dt 


= I + Ah + HAh) 2 +■ 


In all applications discussed in this paper, a discrete form of equation (3) will be used. 
That solution may have some value if an approximation can be found that is efficient, 
accurate and robust. 


Section 2.0. Approximation of the A(t) matrix for an LTV system. Consider the 
discrete analog to the continuous functions above. 


Xn=x n (t m -h) , xl=x n (t m ) , 

® n,l ~ ^0 > ^ n,l ’ 

K) = Kl (tm ~ k ) ’ b n,I = K,l (tm ) , 

= u n(t m -h) , u° n =u n (t m ) , 


xl=x„(t m +h) x[ =x n (t m +rh) 

a \,i = a n,i (t m + h) , a r nJ = a nJ (t m + rh) 
K ,1 = Kl (tm + h) , b r nJ = b n l (t m + rh) 
K=u n (t m +h) u' n = u n (t m + rh) 


The numerical solution of equation (1 ) can be obtained by starting with a term by term 
polynomial approximation of the matrix A(t) 


i(*) = S 


C t,njT 


i - 1 


!h 


i - 1 


r = t-t m , re (0 ,h) 


(4) 
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for schemes of various order. For instance, one can define the 0( A ) scheme where 
A = / - 1 . The coefficients in equation (4) are constructed from discrete values of the 
terms of matrix A(t), that is c inl = c inl (a nl ) . These interpolation polynomials 

can be easily evaluated since Aft) is known for time steps at least up to t m +h. To 
illustrate the use of equation (4 ), the a n l term of the continuous matrix A (t) can be 

approximated by polynomial expressions of orders zero, one and two for re (0 ,h) by the 
following 


0(0) scheme: 
0(1) scheme: 
0(11) scheme: 


a n ,l<j) = C Xn,l where C ln,l =</ 

a n ,i (?) = C, nJ + c 2 nl T / h where c, nJ = a° tJ 

a nj(?) = C ln,l + C 2n,lt / h + C 3 „jT !k 


■'In, l 


= a 


n,l 


-a 


0 

n,l 


where 


■'In, l 


= a 


n,l 


C 2 n,l = 7 ( a h ~ a tl ) ’ C 2n,l ~ 2 (</ " 2a lj + a n,l ) 


In a similar manner, higher order approximations of the matrix Aft) can be constructed 
term by term. The accuracy of this scheme is dependent on the accuracy of the 
approximation of matrix Aft) since equation (3) is exact. Under certain circumstances a 
method based on equation (3) will be an exact numerical scheme. 

Section 2.1. Approximation of the STM for an LTV system. Using the results of 
section 2.0, an algorithm is constructed to compute successive terms of the Peano-Baker 
formula. Since polynomials of order A = / - 1 are used to approximate the matrix A(t) 
term by term, the integration of those expressions over (t m ,t m + h) can be accomplished 
in a straight forward manner. The first two terms are easily found. The computation 
begins as follows: 


(<;)i = c t, n ;/i > i = ll 

, n = \,N , 

1 = 1, N 

(5 a) 

I 

n + h’t m ) = S„,l +^X (<7 ".') 1 

, n = \,N , 

1 = 1, N 

(5 b) 


;=i 


where S n l =1 if n = / and 8 n , = 0 otherwise. The superscript i in equations (5a-b ) is an 

index ranging from 1 to I. In equation (4) the same index i is used as an exponent. 
Equations (5 a) and (5 b) compute the sum of the first two terms in the series equation (3) 
and provide the start of the STM. The terms (<r' , ), computed in equation (5a) are the 
kernel of a recursion that computes the remaining terms of the STM. 

Starting at the double integral (the third term in equation (3)) the following recursive 
algorithm computes the remaining integral terms and adds their contribution to the STM. 
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For j = 0 ,J 


p N 


k=\ i = 1 


(Kh =L2X„, «**'), kp + j+ \) (V/) 

ffy-P+k+jil-V+I T 
,n,i\ u ij )\ 

k=\ i=l 


Kr ,(/ - 1)+2/ ) 2 =it- M+x 


- P + 1 (j + 2) + 1 

(7 > 1 ) 

« ) 2 = t £ c t , B , «* +2 ), /(/> + j + 2) (7 > 1) 




(6 a) 


£=1 i=l 


27 — 1 + 7 ( 7 - 1 ) 


®,/(<.+M,) = ®,;(f.+A,g + A ;+2 S«) 2 , » = 1 = IN (67) 


£=1 


«;)i =«;) 2 


IN 

(6 b) 


(6 c) 

= 1 ,N 


1 = 1, N 

(6 d ) 

1 = 1, N 

(6 c) 


end j 


Successive terms in the integral series are computed starting at j = 0 with the third term 
(second integral). When j = 0, the outer subscript ( ), denotes values computed in 
equations (5). When / > 0, the outer subscript ( ), denotes values from the previous 
iteration of equations (6). There are a total of J + 3 terms of the Neumann series 
computed by this algorithm, combined from both sets of formulae, equations (5) and (6). 
For instance if J = 0, the first 3 terms are computed, namely the first two obtained from 
the former equations and that with j = 0 from the latter recursive formulae. When 
implemented, several indexing details are important. If 1=1, the index p ranging over 1 
to (1-1) is computed once for p = 1, and the sequence with index p over (1-1) to (j+l)(I-l) 
is computed once for p = 0. 

When 1=1, this solution is just the exponential series for a constant matrix (i.e. exp (Ah )) 
or alternately a zero order approximation of a time varying Aft) matrix. Solution of an 
LTV system at this level of approximation corresponds simply to applying the 
exponential series of the steady state STM updated at each time step. Furthermore, when 
the 0(0) scheme is used (7 = 1 ) and 7 = -1 (i.e. only the first two terms computed), this 
algorithm is the explicit backward Euler method. Since this scheme uses as its basis an 
exact solution, it possible with enough terms to compute solutions using large time step 
sizes, as long as the matrix A (t) is sufficiently well behaved so as to be accurately 
approximated by a polynomial expansion. In the cases computed so far, relatively few 
terms in the series have been required to achieve very good accuracy. 

It would be desirable to obtain a bound on the size of the leading order truncated term 
following the last of the 7 + 3 terms used in the computation. Denote the size of the 
leading order truncated term by Ej+3. While it may be possible to obtain Ej+3 for any time 
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varying system it is likely to be unduly restrictive if based on the most general form of 
the equations above. Rather consider the Lxn order linear time invariant ( LTI ) equation 


d n x l 

df 


Aj -t; i 


/ = U, L 


where A, , l = are the coefficients of a diagonalized system. It has been found 
empirically by the author that the largest value of the highest order term in the series 
truncation error, using J + 3 terms of the series, will have for bound 

7 ./+ 2/0 \<2 T ■ 1 

E J+ 3 < for a = {1 + int( )} (7) 

( J + 2)i n 

where /l max = max^ |,...,|/l z |] . In the examples that follow, the largest magnitude entry 
of the matrix A (t) at each time step will be used as an estimate of /L max . By doing this, 
the assumption is that this will result in a sufficiently accurate estimate of the truncation 
error. Equation (7) does not of course represent the total error in the scheme, which 
would also include the error in the approximation in equation (4). 

The solution of the last term in equation (2) is facilitated by the fact that values of the 
STM at successive time steps are already available, which can be used to numerically 
evaluate the integral term. For constant step size h, write 

£2 =®(t m +hJ m )B(tJu(t m ) 

& = + h, )B(t m _ x )u(t m _ x ) 

= 0(t m +h,t m )Ci 

A recursion for computing the vectors £ l m ,...,^” /+/ can easily be consfructed from these 
relations. The vector of integrals is evaluated by defining I e R v where 

t m +h 

I = j o (t m + h,T)B(T)u(T)dT 


or an approximation of this integral 

0 *=1 


( 8 ) 
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where 1-1 is the order, g = T - t m and the vector c i has the functional dependence 
c t = c i (f] n ,...,^~ m I+2 ) , £,...,£; /+2 e R ,c t e v . To illustrate the procedure, the 
approximation (equation (8)) can be evaluated for orders zero and one by 

0(0) scheme: I ~ hc: x where c x = ° 

0(1) scheme: l~h(c x +\c 2 ) where c, =C , c 2 = £-#° 

This completes development of the method for an LTV system. The next section outlines 
modification of the procedure for an NLTV system. 

Section 2.3. Approximation of the STM for an NLTV system. Consider now the case 
in which 


x = A(t, x)x + B(t)u(t ) (9) 

Suppose that the linear and nonlinear parts of the matrix A( t,x ) can be expressed as 

A(t,x) = A(t) + F(x) 


where 

= a nJ dm ~ h ) + fnj , «nj = « nj (*m ) + fnj 

where 

fn1=fn,Mt m ~h)) , .C =/,,(.V(/„.)) 


Ki = a n^ K +rh) + f nl 
fnj =fnj(x(t m +rh )) 


and the coefficients a\,, ah,... are the time varying part. A local approximation of 

A(t,x ) can be found using equation (4) where now c inl = c inl (a^ ) ■ Because 
the construction of these coefficients is a crucial step additional explanation is required. 
The nonlinear terms a nJ are again approximated by polynomial expressions. For 

example, orders zero and one are given by 


0(0) scheme: 
0(1) scheme: 
where 


a n i(r) = c \ nl 


a n j(T) = c Xn , +c 2 nJ r/h 


where c ln ,=a° n , +/° 


C 1 n,l u n,l T J nJ ’ c 2 n,l u n,l u n,l T J n,l J nl 
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and re (0 ,h) . The value of the solution x(t m + h) has not yet been computed and 
therefore the function /' ' ; is not yet known. This tenn is computed by extrapolation, that 
is f' n [ = f n d (x 1 ), where x 1 = g(x° ,...,x~ I+l ) is derived from previous solutions. The 
rest of the method follows that outlined in the previous section. 


Section 3.0. Examples. In this section solutions of LTV and NLTV systems are 
computed. 

Example 3.1. This example computes the following second order system 

x = t 4 x , x(0) = 0, x(0) = 1 
or 


where 


A = 


0 1 

t 4 0 


computed on t : 0 — > 1 . The exact solution is given by the series 


x(t) = x n 


v*-i (6m - 5 )( 6m - 6 ) 

1 


.6k - 5 


n 


k=\ \m- 


= i (6m -5)(6m -6) 


{6k - 5) t 


6k-6 


Figure 1 presents the error in the 0(1), 0(11), O(III) and O(IV) schemes compared with 
the second order implicit Euler method. Convergence to the exact solution is shown as 
the order of the approximation of A(t) increases. The O(IV) scheme is exact to machine 
accuracy since the time variation of the matrix A(t) is exactly captured with a 4 lh order 
polynomial. 

This example illustrates one of the advantages of this method. If the time variation is 
sufficiently modeled with the order polynomial used, the method is nearly insensitive to 
the size of the time step. Using the 0( IV) scheme in this example, any time step size is 
permissible since the exact solution is being computed. This aspect will be explored more 
fully in the next example. 
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q 2nd ord. implicit Euler, h=0.04 

A 0(1) scheme, h=0.04 

9 0(11) scheme, h=0.04 

q 0(111) scheme, h=0.04 



Figure 1. Percent error. 


Example 3.2. This example computes the response of a pitch/plunge apparatus initially 
oscillating as a sinusoid, due to a change in the torsion stiffness to % its original value. 
The plunge and initial torsion stiffness are K h =1.21x10 6 N / m and K aX =6.68x10’ 

N -ml rad . The mass is m h = 26.64 kg , the mass pitch moment of inertia is 
I c/A = 0.086 kg - m 2 , and the static offset is s a = 0.378 kg-m. These values are for an 
actual system set up to study transonic flutter. 3 The functional variation of the torsion 
stiffness is given by 


K a =K a] (\-f) + K a2 f tl <t<t l+ At 

where 

/ = y(l-cos[;r(/-t 1 )/Ai]) 

and 6 = 0.24 sec, At = 0.01 sec. The natural frequencies of the plunge and torsion 
modes are 205.4 rad/sec and 299.3 rad/ sec for t < i, . The second order system is 
diagonalized, resulting in the equation 


q = Aq + Bu . 


The matrix A(t) has the values 
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-42189 

0 



"-39298 

10533" 



for t <t x . 

A = 



0 

-89580 



10533 

-50876_ 


for t>r l + A t 


and 


B = 


1 0 
0 1 


u = 


The generalized variables q x and q 2 roughly correspond to pitch and plunge. These 
equations are written in state variable form and solved using equations (5) and (6) with a 
sufficient number of terms adjusted at each time step to ensure E J+3 <1x1 (T 7 (see 
equation (7)). 

Computations are performed with h = 0.0005, 0.001, 0.002, and 0.004 sec (approximately 
64 to 8 time steps/pitch cycle). The ensuing results are compared with an integration 
using a variable time step size 4 th order Runge-Kutta (R-K) scheme with tolerance of 

I x 10“ 7 . The resulting time traces are shown in Figures 2-1 1 . Figures 2-7 are computed 
with zero input, i.e. u x =u 2 = 0 , but with initial conditions q A (0) = q 2 (0) = 1 . Figures 8- 

I I are computed with initial conditions q x (0) = 0, q 2 (0) = 0 , but with an impulse input 
for u 2 centered at t = . 1 sec . 

Figure 2 gives an overview of the time history of variable q 2 solved with the R-K 
method, showing the torsion stiffness change at time /, = 0.24 sec and continuing 
afterwards for about 40 cycles. The R-K solution required an average time step size of 
have = 0-0001 sec to achieve the desired error tolerance and about 1 0000 time steps over 
the interval shown. 

Figure 3 compares the O(III) scheme solution at h = 0.004 sec with the R-K simulation. 
Even after nearly 50 cycles of oscillation there is no discernable difference in the 
solutions. A solution computed with the 2 nd order implicit Euler method using 
h = 0.002 sec is compared with the 4 th order variable time step R-K solution in Figure 4. 
Although the very large phase error of the Euler method for this problem makes that 
solution unusable, it does illustrate the difficulty in obtaining accurate solutions over a 
very large time span. By contrast, all of the solutions shown using the present scheme are 
quite accurate. 

The solutions in Figures 5-7 using the 0(1) to the O(III) scheme are computed at 
successively smaller time step sizes. Figure 5 presents solutions using the 0(1) scheme at 
successive time step sizes. Figure 6 presents solutions using the O(II) scheme, while 
Figure 7 shows solutions using the O(III) scheme, each at successively smaller time step 
sizes. Each set of solutions demonstrates the convergence of the present schemes to the 
4 th order R-K solution as the time step size is decreased. The 0(1) solution in Figure 5 
converges, albeit slowly. As expected, the solution using the 0( III) scheme seen in 
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Figure 2. Second mode response to smooth change in torsion stiffness. 

Figure 7, converges much more rapidly with time step reduction to the R-K solution than 
the solution using the 0(1) scheme seen in Figure 5. 

The last simulations using the pitch/plunge set up, seen in Figures 8-11, show the 
response due to an impulse input. Figure 8 shows the overall solution and the impulse 
input that initiates the dynamic response. The discrete impulse input centered at 
t = 0.1 sec is shown in Figure 9 at successive time step sizes. That figure illustrates the 
coarseness of the approximation of the impulse at the largest time step. Even at the 
largest time step size, at which the impulse is defined by only three time steps, the 
relative insensitivity of the O(III) scheme to time step size can be seen in the results of 
Figures 10-11. 

The last three examples apply the method to NLTV systems. 

Example 3.3. The equation 

x, = x 2 , x 2 = -cx, + g cos ht - dx 2 - ax" (10) 

is solved. When n = 2 this equation is the Flelmholtz oscillator arising in the modeling of 
ship capsizing. 4 When n = 3 the Duffing oscillator can be modeled, arising in the study of 
electronic oscillators 5 , or structural dynamics in which there is nonlinear stiffness. 6 The 
parameters are c = -1, a = b = 1, g = 0.3, d = 0.15, n = 3 . 
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2 r 


4th order variable t.s. R-K (MATLAB) 
0(111) scheme, h = 0.004, 8 1. s./cycle 


□ 



Figure 3. Second mode response to smooth change in torsion stiffness. 

4th order variable t.s. R-K (MATLAB) 



Figure 4. Second mode response to smooth change in torsion stiffness. 
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Figure 5. Second mode response to smooth change in torsion stiffness. 



Figure 6. Second mode response to smooth change in torsion stiffness. 
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Figure 7. Second mode response to smooth change in torsion stiffness. 



Figure 8. Impulsive input and second mode response to input U 2 (t). 
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Figure 1 1. Second mode response to input U 2 (t). 


Solutions of this equation using the 0(0) - 0( IV) schemes with initial condition 

x,(0) = -1, x 2 (0) = 1 are in Figure 12 a)-e). The tolerance in the present schemes is set at 

E J+ 3 < lxlO 7 . A solution is also computed with the variable order multi-step solver 

odel5s in MATLAB, with tolerance of 1 x 10" 7 . The variable order solution takes about 
4000 time steps over the interval shown, with an average time step size of 0.012. As the 
order of the present scheme is increased, the solution converges to the reference solution: 
as shown in Figure 12 e) the O(TV) scheme with h = 0.05 is identical to the reference 
solution. Flowever, as shown in Figure 13, with just a slight change in the initial 
condition to x, (0) = -1 , x 2 (0) = 1.001 the O(IV) scheme with h = 0.05 shows divergence 
from the reference solution toward the end of the simulation. Time step reduction to h = 
0.01 is required using the O(IV) scheme for complete convergence. 

Example 3.4. Equation (10) is solved with values now given by c = -y , a = j , 

b = 0.79 , g = 0.095 , d = 0.1 , n = 3 . This is the twin well oscillator with chaotic motion 

between attractors. 5 For the reference solution, the MATLAB ode 1 5s solver is employed. 

Since this case is much more sensitive to numerical error than the last, the tolerance is 

now set at lxl0~ 14 . At this level of tolerance, the average time step size is 

h me = 0.0005 with 600,000 time steps required. This solution and that computed using 

the O(V) scheme with h = 0.002, and E J+3 < lxlO -12 are shown in Figure 14. The O(V) 
scheme and the MATLAB solution diverge only slightly and only toward the end of the 
simulation. 
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Figure 12 a)-e). Numerical solution of the Duffing equation, x(0)=[-l ,1 .]. 
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MATLAB variable order solver (ode15s) 



Figure 13. Numerical solution of the Duffing equation, x(0)=[-l, 1.001]. 



Figure 14. Numerical solution of the twin well equation, x(0) = [- 1 ,-0.5] . 
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Example 3.5. The van der Pol equation 


X = //( I - x 2 )x - X 

is computed with // = 10. At this value of p , the equation is only moderately stiff. 
Figure 15 shows the solution in the region t = 6-11 in which a transition in the solution 
occurs; before t = 6 the solutions are nearly constant and are identical. The result 
obtained from the O(III) scheme, with h = 0.05 and A /+3 < lxlO -6 , is compared with that 
from a variable order multi-step solver (odel5s in MATLAB) with tolerance set at 
1x10'°. Slight differences in the solutions appear only in the area in which there are 
large solution gradients. The nearly identical results indicate that, at the large time step 
used, the present O(III) scheme computes the rapidly changing solution well. 



Figure 15. Numerical solution of the van der Pol equation, p = 10. 

Section 4.0. Conclusions. A variable order method of integrating initial value ordinary 
differential equations that is based on the state transition matrix (STM) has been 
developed. Although the STM solution formally is valid only for linear equations, the 
discrete numerical method can be applied to both linear time variant and to nonlinear 
systems of equations. When the time variation of a linear system can be exactly modeled 
by a polynomial expansion, the scheme will give an exact result for large time step size. 
Accurate solutions in all the cases shown can be achieved for linear time variant systems 
at large time step sizes. The type of nonlinear system is limited to those that can be 
modeled within a state space frame work, that is, in which the nonlinearity is contained in 
the A(t,x), B(t,x) matrix terms. Several nonlinear problems such as the Duffing equation 
and twin well attractor, exhibiting chaotic behavior, have also been computed. Higher 
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order solutions of these equations compare well with reference solutions performed with 
the standard 4 th order Runge-Kutta and variable order methods. 

The present method has been developed for integration of modal structural, flight 
dynamics and control system equations in conjunction with a computational fluid 
dynamics solver, an application for which it is ideally suited. Applications might include 
the simulation of structural dynamics with time varying parameters such as stiffness, 
mass or geometry, simulation of a fully adaptive control system, or flight dynamics with 
variation in flight parameters. Other applications are also being considered. Apart from 
its use in a numerical integration scheme, it offers a generalized numerical technique for 
the computation of the STM for any time varying linear system, among other possible 
applications. 
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